{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n",
      "15\n",
      "16\n",
      "17\n",
      "18\n",
      "19\n",
      "20\n",
      "21\n",
      "22\n",
      "23\n",
      "24\n",
      "25\n",
      "26\n",
      "27\n",
      "28\n",
      "29\n",
      "30\n",
      "31\n",
      "32\n",
      "33\n",
      "34\n",
      "35\n",
      "36\n",
      "37\n",
      "38\n",
      "39\n",
      "40\n",
      "41\n",
      "42\n",
      "43\n",
      "44\n",
      "45\n",
      "46\n",
      "47\n",
      "48\n",
      "49\n",
      "50\n",
      "51\n",
      "52\n",
      "53\n",
      "54\n",
      "55\n",
      "56\n",
      "57\n",
      "58\n",
      "59\n",
      "60\n",
      "61\n",
      "62\n",
      "63\n",
      "64\n",
      "65\n",
      "66\n",
      "67\n",
      "68\n",
      "69\n",
      "70\n",
      "71\n",
      "72\n",
      "73\n",
      "74\n",
      "75\n",
      "76\n",
      "77\n",
      "78\n",
      "79\n",
      "80\n",
      "81\n",
      "82\n",
      "83\n",
      "84\n",
      "85\n",
      "86\n",
      "87\n",
      "88\n",
      "89\n",
      "90\n",
      "91\n",
      "92\n",
      "93\n",
      "94\n",
      "95\n",
      "96\n",
      "97\n",
      "98\n",
      "99\n",
      "100\n",
      "101\n",
      "102\n",
      "103\n",
      "104\n",
      "105\n",
      "106\n",
      "107\n",
      "108\n",
      "109\n",
      "110\n",
      "111\n",
      "112\n",
      "113\n",
      "114\n",
      "115\n",
      "116\n",
      "117\n",
      "118\n",
      "119\n",
      "120\n",
      "121\n",
      "122\n",
      "123\n",
      "124\n",
      "125\n",
      "126\n",
      "127\n",
      "128\n",
      "129\n",
      "130\n",
      "131\n",
      "132\n",
      "133\n",
      "134\n",
      "135\n",
      "136\n",
      "137\n",
      "138\n",
      "139\n",
      "140\n",
      "141\n",
      "142\n",
      "143\n",
      "144\n",
      "145\n",
      "146\n",
      "147\n",
      "148\n",
      "149\n",
      "150\n",
      "151\n",
      "152\n",
      "153\n",
      "154\n",
      "155\n",
      "156\n",
      "157\n",
      "158\n",
      "       Beijing_S Beijing_E  Beijing_I   Beijing_R Beijing_D   Beijing_H  \\\n",
      "0    2.15463e+07     0.059  0.0865333           0         0           0   \n",
      "1    2.15506e+07   0.13069   0.176054  0.00454243         0  0.00378205   \n",
      "2    2.15548e+07  0.215445   0.269724   0.0137436         0   0.0121596   \n",
      "3    2.15591e+07  0.313748   0.368697   0.0277614         0   0.0259702   \n",
      "4    2.15634e+07   0.42619   0.474134   0.0467952         0   0.0460822   \n",
      "..           ...       ...        ...         ...       ...         ...   \n",
      "154  2.17244e+07   13.5006    63.7135     284.522         0     331.183   \n",
      "155  2.17244e+07   13.3198    62.9066     286.194         0     332.049   \n",
      "156  2.17244e+07   13.1418    62.1094     287.845         0     332.903   \n",
      "157  2.17244e+07   12.9666    61.3216     289.475         0     333.745   \n",
      "158  2.17244e+07   12.7941    60.5433     291.085         0     334.576   \n",
      "\n",
      "         Wuhan_S  Wuhan_E  Wuhan_I  Wuhan_R Wuhan_D  Wuhan_H  \n",
      "0    1.08886e+07  168.445  223.755  5.77428       1  50.6154  \n",
      "1    1.08843e+07  187.387  228.591  11.6448       1  61.4131  \n",
      "2      1.088e+07  206.914  234.514    17.64       1  73.4251  \n",
      "3    1.08756e+07  227.116   241.53  23.7883       1  86.6888  \n",
      "4    1.08713e+07  248.084  249.654  30.1183       1  101.248  \n",
      "..           ...      ...      ...      ...     ...      ...  \n",
      "154   1.0672e+07  1035.39  6569.51  30937.6       1  37345.1  \n",
      "155   1.0672e+07  1016.81  6463.45    31110       1  37411.5  \n",
      "156  1.06719e+07  998.645  6358.98  31279.6       1  37476.7  \n",
      "157  1.06719e+07  980.872   6256.1  31446.5       1  37540.7  \n",
      "158  1.06718e+07   963.48  6154.77  31610.7       1  37603.6  \n",
      "\n",
      "[159 rows x 12 columns]\n",
      "Running time: 1466.156789779663 Seconds\n"
     ]
    }
   ],
   "source": [
    "#本代码作者 宋鹏飞 西安交通大学 数学与统计学院\n",
    "#时间2020年03月08日\n",
    "#本代码为突发疾病多阶段多尺度模型核心算法， 以武汉市和北京市为例\n",
    "#由于作者没有受过系统的编程训练，纯属业余，代码过于粗糙，有编程‘洁癖’的敬请谅解\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import time\n",
    "start=time.time() #计算模拟时间 \n",
    "#初始时间1月9日24时\n",
    "Tspan=53 #初始时刻, 可修改\n",
    "dt=1/3 #单位时间, 可修改\n",
    "Tlock=14 #封城时间 1月23日封城\n",
    "\n",
    "Tbeijinghome=15 #1月24日一级响应\n",
    "Tbeijingrelease=32  #逐渐恢复市内交通\n",
    "Twuhanhome=14  #武汉home\n",
    "Thubeirelease=114 #湖北五月初恢复交通\n",
    "\n",
    "Trelease=145 #湖北六月初离汉通道开启\n",
    "#城市属性  可自行修改城市及区\n",
    "Citys=['Beijing','Wuhan']\n",
    "Wuhanqu=['江岸区','江汉区','硚口区','汉阳区','武昌区','青山区',\n",
    "         '洪山区','蔡甸区','江夏区','黄陂区','新洲区','东西湖区','汉南区']\n",
    "Beijingqu=['密云区','延庆区','朝阳区','丰台区','石景山区','海淀区','门头沟区', \n",
    "           '房山区','通州区','顺义区','昌平区','大兴区','怀柔区','平谷区','东城区','西城区']\n",
    "\n",
    "Wuhan_renkou=10000*np.array([96.24,72.96,86.85,65.27,127.63,52.88,163.75,\n",
    "                       73.50,91.37,98.83,90.21,56.25,13.55]) #武汉人口\n",
    "\n",
    "WUHAN=pd.DataFrame(data=Wuhan_renkou,index=Wuhanqu,columns=['Renkou'])\n",
    "Beijing_renkou=10000*np.array([49.5,34.8,360.5,210.5,59.0,335.8,33.1,118.8,157.8,116.9,210.8,179.6,\n",
    "                        41.4,45.6,82.2,117.9]) #北京人口\n",
    "BEIJING=pd.DataFrame(data=Beijing_renkou,index=Beijingqu,columns=['Renkou'])\n",
    "Patches={'Beijing':Beijingqu,'Wuhan':Wuhanqu}\n",
    "\n",
    "#Wuhancases=np.array([258,363,425,495,572,618,698,1590,1905,2261,2639,3215,4109,5142,6384,8351,10117,11618,13603,14982,16902,18454,19558,...\n",
    "   #20630+12364,21960+14031,37914,39462,41152,42752,44412,45027,45346,45660])\n",
    "\n",
    "Atr=pd.DataFrame(index=Citys)#城市属性\n",
    "Atr['Pnum']=np.zeros(len(Citys))\n",
    "for i in Citys:#添加城市属性区的数目\n",
    "    Atr.loc[i,'Pnum']=len(Patches[i])\n",
    "Atr['Province']=['Beijing','Hubei'] #添加城市省份\n",
    "\n",
    "\n",
    "# 按照模型将人群分为c*w*h个子类，c所在城市, w工作所在地， h为家庭所在地\n",
    "INDEX=[]\n",
    "for c in Citys:\n",
    "    for w in Patches[c]:\n",
    "        for h in Patches[c]:\n",
    "            INDEX.append(c+h+w)\n",
    "            \n",
    "            \n",
    "# 输入的每个人群子类的参数属性\n",
    "PARA=pd.DataFrame(index=INDEX,columns=['Cityname','S_initial','E_initial',\n",
    "                           'I_initial','R_initial','D_initial','H_initial'])\n",
    "for c in Citys:#更新城市名字\n",
    "    for i in PARA.index:\n",
    "        if c in i:\n",
    "            PARA.loc[i,'Cityname']=c\n",
    "PARA['S_initial']=np.ones(len(INDEX))\n",
    "c='Wuhan'\n",
    "for h in Patches[c]:\n",
    "    for w in Patches[c]:\n",
    "        if h==w:\n",
    "            PARA.loc[c+h+w,'S_initial']=WUHAN.loc[h,'Renkou']/3+WUHAN.loc[h,'Renkou']**2/sum(Wuhan_renkou)*2/3\n",
    "        else:\n",
    "            PARA.loc[c+h+w,'S_initial']=2/3*WUHAN.loc[h,'Renkou']*WUHAN.loc[w,'Renkou']/sum(Wuhan_renkou)\n",
    "c='Beijing'\n",
    "for h in Patches[c]:\n",
    "    for w in Patches[c]:\n",
    "        if h==w:\n",
    "            PARA.loc[c+h+w,'S_initial']=BEIJING.loc[h,'Renkou']/3+BEIJING.loc[h,'Renkou']**2/sum(Beijing_renkou)*2/3\n",
    "        else:\n",
    "            PARA.loc[c+h+w,'S_initial']=2/3*BEIJING.loc[h,'Renkou']*BEIJING.loc[w,'Renkou']/sum(Beijing_renkou)    \n",
    "'''AAA=0\n",
    "c='Wuhan'\n",
    "for h in Patches[c]:\n",
    "    for w in Patches[c]:\n",
    "        AAA=AAA+PARA.loc[c+h+w,'S_initial']\n",
    "print(AAA)'''\n",
    "PARA['E_initial']=np.zeros(len(INDEX))\n",
    "PARA.loc['Wuhan江汉区江汉区','E_initial']=150\n",
    "PARA['I_initial']=np.zeros(len(INDEX))\n",
    "PARA.loc['Wuhan江汉区江汉区','I_initial']=220\n",
    "PARA['R_initial']=np.zeros(len(INDEX))\n",
    "PARA['D_initial']=np.zeros(len(INDEX))\n",
    "PARA.loc['Wuhan江汉区江汉区','D_initial']=1\n",
    "PARA['H_initial']=np.zeros(len(INDEX))\n",
    "PARA.loc['Wuhan江汉区江汉区','H_initial']=41\n",
    "#print(PARA)\n",
    "\n",
    "\n",
    "# 输出\n",
    "#各个地方的S， E, I, R^{c,w,h}\n",
    "COLUMNS=[] #每列\n",
    "for i in INDEX:\n",
    "    A=[i+'_S',i+'_E',i+'_I',i+'_R',i+'_D',i+'_H']\n",
    "    COLUMNS.extend(A)\n",
    "#print(COLUMNS)\n",
    "Solution=pd.DataFrame(columns=COLUMNS)\n",
    "\n",
    "\n",
    "#定义控制参数C\n",
    "def Control(c,i,t): #city, location, time, initial time is Jan 9th,24:00 OK\n",
    "    cnh=20#控制措施产生效果\n",
    "    cnhk=0.14#效果\n",
    "    cnhday=5#花了多少天\n",
    "    cwh=20#控制措施产生效果\n",
    "    cwhday=8#花了多少天\n",
    "    cwhk=0.083#效果\n",
    "    def cNH(t):#湖北以外控制\n",
    "        if t<=cnh: #1月30 号之后\n",
    "            c=1\n",
    "        elif t>cnh+cnhday:\n",
    "            c=cnhk  #可调\n",
    "        else:\n",
    "            c=np.exp(np.log(cnhk)*(t-cnh)/cnhday)\n",
    "        return c\n",
    "    def cH(t):#湖北控制\n",
    "        if t<=cwh: #2月3号之后\n",
    "            c=1\n",
    "        elif t>cwh+cwhday:\n",
    "            c=cwhk  #可调\n",
    "        else:\n",
    "            c=np.exp(np.log(cwhk)*(t-cwh)/cwhday)\n",
    "        return c\n",
    "    if Atr.loc[c,'Province']=='Hubei':\n",
    "        control=cH(t)\n",
    "    else:\n",
    "        control=cNH(t)\n",
    "    return control\n",
    "#print(Control('Wuhan',16,26))\n",
    "#定义潜伏感染者感染率\n",
    "def Betae(c,i,t):#Exposed 感染率， i表示location, t时间 OK\n",
    "    k=0.95\n",
    "    if Atr.loc[c,'Province']=='Hubei':\n",
    "        betae=0.24*k\n",
    "    else:\n",
    "        betae=0.20*k\n",
    "    return betae\n",
    "#print(Betae('Beijing',10))\n",
    "#定义有症状感染者感染率\n",
    "def Betai(c,i,t):#Infected 感染率， i表示location, t时间 OK\n",
    "    k=0.95\n",
    "    if Atr.loc[c,'Province']=='Hubei':\n",
    "        betai=0.24*k\n",
    "    else:\n",
    "        betai=0.20*k\n",
    "    return betai\n",
    "#定义市内出行概率\n",
    "def Prob(c,w,h,i,t): # c city, w work place, h, home place, i, location, t time OK\n",
    "    if c=='Beijing':\n",
    "        if t<=Tbeijinghome:   \n",
    "            prob=1/Atr.loc[c,'Pnum']\n",
    "        elif t>Tbeijinghome and t<Tbeijingrelease:\n",
    "            if i==h:\n",
    "                prob=1\n",
    "            else:\n",
    "                prob=0\n",
    "        else:\n",
    "            if (np.floor(t)-Tbeijingrelease+1)%7==5 or (np.floor(t)-Tbeijingrelease+1)%7==6:#双休日\n",
    "                if t-np.floor(t)>0 and t-np.floor(t)<=1/3+0.001:#白天\n",
    "                    prob=1/Atr.loc[c,'Pnum']\n",
    "                else:\n",
    "                    if i==h:\n",
    "                        prob=1\n",
    "                    else:\n",
    "                        prob=0\n",
    "            else:\n",
    "                if t-np.floor(t)>0 and t-np.floor(t)<=1/3+0.001:#白天\n",
    "                    if i==w:\n",
    "                        prob=1\n",
    "                    else:\n",
    "                        prob=0\n",
    "                else:\n",
    "                    if i==h:\n",
    "                        prob=1\n",
    "                    else:\n",
    "                        prob=0\n",
    "    if c=='Wuhan':\n",
    "        if t<=Twuhanhome:   \n",
    "            prob=1/Atr.loc[c,'Pnum']\n",
    "        elif t>Twuhanhome and t<Thubeirelease:\n",
    "            if i==h:\n",
    "                prob=1\n",
    "            else:\n",
    "                prob=0\n",
    "        else:\n",
    "            if (np.floor(t)-Thubeirelease+1)%7==5 or (np.floor(t)-Thubeirelease+1)%7==6:#双休日\n",
    "                if t-np.floor(t)>0 and t-np.floor(t)<=1/3+0.001:#白天\n",
    "                    prob=1/Atr.loc[c,'Pnum']\n",
    "                else:\n",
    "                    if i==h:\n",
    "                        prob=1\n",
    "                    else:\n",
    "                        prob=0\n",
    "            else:#工作日\n",
    "                if t-np.floor(t)>0 and t-np.floor(t)<=1/3+0.001:#白天\n",
    "                    if i==w:\n",
    "                        prob=1\n",
    "                    else:\n",
    "                        prob=0\n",
    "                else:\n",
    "                    if i==h:\n",
    "                        prob=1\n",
    "                    else:\n",
    "                        prob=0            \n",
    "    return prob\n",
    "'''for i in range(180):\n",
    "    print(Prob('Wuhan','汉南区','江汉区','汉南区',i/3))'''\n",
    "#定义城市与城市之间迁徙\n",
    "def Diffusion(a,b,t): # a city to b city, distance induced diffusion rate OK\n",
    "    if t>Tlock and t<Trelease:\n",
    "        diffusion=0\n",
    "    else:\n",
    "        diffusion=0.00118\n",
    "    return diffusion\n",
    "def Migration(a,b,t):#a city to b city,migration probability OK\n",
    "    if a=='Wuhan' and b=='Beijing':\n",
    "        migration=1\n",
    "    else:\n",
    "        migration=0\n",
    "    return migration\n",
    "#定义迁入人群概率Q，减模型\n",
    "def Qprob(c,w,h,t):# c city, w work place, h, home place, t time OK\n",
    "    qprob=1/Atr.loc[c,'Pnum']**2\n",
    "    return qprob\n",
    "#定义恢复率\n",
    "def Gamma(c,i,t): # c city, i patch, t time OK\n",
    "    y=1/12.7\n",
    "    return y\n",
    "#定义因病死亡率\n",
    "def Alpha(c,i,t): # c city, i patch, t time OK\n",
    "    a=0;\n",
    "    return a\n",
    "#定义潜伏期长度\n",
    "def Sigma(): #OK\n",
    "    s=1/5.2\n",
    "    return s\n",
    "\n",
    "#定义中间变量\n",
    "COLUMNS1=[] #每列\n",
    "for i in INDEX:\n",
    "    A=[i+'_S1',i+'_E1',i+'_I1',i+'_R1',i+'_D1',i+'_H1']\n",
    "    COLUMNS1.extend(A)\n",
    "COLUMNS2=[] #每列\n",
    "for i in INDEX:\n",
    "    A=[i+'_S2',i+'_E2',i+'_I2',i+'_R2',i+'_D2',i+'_H2']\n",
    "    COLUMNS2.extend(A)\n",
    "#print(COLUMNS1)\n",
    "Midpar1=pd.DataFrame(data=np.zeros((1,len(COLUMNS1))),columns=COLUMNS1)\n",
    "Midpar2=pd.DataFrame(data=np.zeros((1,len(COLUMNS2))),columns=COLUMNS2)\n",
    "\n",
    "#print(type(Midpar2.iloc[:,0]))\n",
    "#代入初值\n",
    "#!!!!!!!算法开始，微分方程采取最简单的差分，可自行修改微分方程解法\n",
    "for i in INDEX:\n",
    "    for j in ['S','E','I','R','D','H']:\n",
    "        Midpar1[i+'_'+j+'1']=PARA.loc[i,j+'_initial']\n",
    "#print(Midpar1)\n",
    "for i in range(int(Tspan/dt)):\n",
    "    print(i)\n",
    "    for c in Citys:\n",
    "        #计算S_{i}^{c} E_{i}^{c} I R_{i}^{c}\n",
    "        Nlocation=pd.DataFrame(data=np.zeros((len(Patches[c]),5)),\n",
    "                          columns=['Ssum','Esum','Isum','Rsum','Nsum'],index=Patches[c])\n",
    "        #print(type(Nlocation.loc[patch,'Ssum']))\n",
    "        for patch in Patches[c]:\n",
    "            for hh in Patches[c]:\n",
    "                for ww in Patches[c]:\n",
    "                    Suma=Nlocation.loc[patch,'Ssum']+Prob(c,ww,hh,patch,i*dt)*Midpar1.loc[0,c+hh+ww+'_S1']\n",
    "                    #print(Suma)\n",
    "                    #print(type(Prob(c,ww,hh,patch,i*dt)*Midpar1.loc[0,c+hh+ww+'_S1']))\n",
    "                    Sumb=Nlocation.loc[patch,'Esum']+Prob(c,ww,hh,patch,i*dt)*Midpar1.loc[0,c+hh+ww+'_E1']\n",
    "                    Sumc=Nlocation.loc[patch,'Isum']+Prob(c,ww,hh,patch,i*dt)*Midpar1.loc[0,c+hh+ww+'_I1']\n",
    "                    Sumd=Nlocation.loc[patch,'Rsum']+Prob(c,ww,hh,patch,i*dt)*Midpar1.loc[0,c+hh+ww+'_R1']\n",
    "                    Nlocation.loc[patch,'Ssum']=Suma\n",
    "                    Nlocation.loc[patch,'Esum']=Sumb\n",
    "                    Nlocation.loc[patch,'Isum']=Sumc\n",
    "                    Nlocation.loc[patch,'Rsum']=Sumd\n",
    "                    Nlocation.loc[patch,'Nsum']=Suma+Sumb+Sumc+Sumd\n",
    "        Ntotal=sum(Nlocation['Nsum'].values)\n",
    "        #print(Ntotal) OK\n",
    "        #计算扩散项:迁出 和迁入 \n",
    "        imigration_rate=0\n",
    "        emigration_num_S=0\n",
    "        emigration_num_E=0\n",
    "        emigration_num_I=0\n",
    "        emigration_num_R=0\n",
    "        for b in Citys:\n",
    "            imigration_rate=imigration_rate+Migration(c,b,i*dt)*Diffusion(c,b,i*dt)\n",
    "            #计算 S^{b}\n",
    "            SSUM=0\n",
    "            ESUM=0\n",
    "            ISUM=0\n",
    "            RSUM=0\n",
    "            for hb in Patches[b]:\n",
    "                for wb in Patches[b]:\n",
    "                    SSUM=SSUM+Midpar1.loc[0,b+hb+wb+'_S1']\n",
    "                    #print(SSUM)\n",
    "                    ESUM=ESUM+Midpar1.loc[0,b+hb+wb+'_E1']\n",
    "                    ISUM=ISUM+Midpar1.loc[0,b+hb+wb+'_I1']\n",
    "                    RSUM=RSUM+Midpar1.loc[0,b+hb+wb+'_R1']\n",
    "            emigration_num_S=emigration_num_S+Migration(b,c,i*dt)*Diffusion(b,c,i*dt)*SSUM\n",
    "            emigration_num_E=emigration_num_E+Migration(b,c,i*dt)*Diffusion(b,c,i*dt)*ESUM\n",
    "            emigration_num_I=emigration_num_I+Migration(b,c,i*dt)*Diffusion(b,c,i*dt)*ISUM\n",
    "            emigration_num_R=emigration_num_R+Migration(b,c,i*dt)*Diffusion(b,c,i*dt)*RSUM\n",
    "        for h in Patches[c]:\n",
    "            for w in Patches[c]:       \n",
    "                Infection=0\n",
    "                Recover=0\n",
    "                Dead=0\n",
    "                for patch in Patches[c]:\n",
    "                    I_sum=Nlocation.loc[patch,'Isum']\n",
    "                    E_sum=Nlocation.loc[patch,'Esum']\n",
    "                    N_sum=Nlocation.loc[patch,'Nsum']\n",
    "                    Infection=Infection+Prob(c,w,h,patch,i*dt)*Control(c,patch,i*dt)*(Betai(c,patch,i*dt)*I_sum/N_sum+Betae(c,patch,i*dt)*E_sum/N_sum)\n",
    "                    Recover=Recover+Prob(c,w,h,patch,i*dt)*Gamma(c,patch,i*dt)\n",
    "                    Dead=Dead+Prob(c,w,h,patch,i*dt)*Alpha(c,patch,i*dt)\n",
    "                #print(Recover) OK\n",
    "                qprob=(Midpar1.loc[0,c+h+w+'_S1']+Midpar1.loc[0,c+h+w+'_E1']+Midpar1.loc[0,c+h+w+'_I1']+Midpar1.loc[0,c+h+w+'_R1'])/Ntotal\n",
    "                Midpar2.loc[0,c+h+w+'_S2']=Midpar1.loc[0,c+h+w+'_S1']+dt*(-Midpar1.loc[0,c+h+w+'_S1']*Infection-imigration_rate*Midpar1.loc[0,c+h+w+'_S1']+qprob*emigration_num_S)\n",
    "                Midpar2.loc[0,c+h+w+'_E2']=Midpar1.loc[0,c+h+w+'_E1']+dt*(Midpar1.loc[0,c+h+w+'_S1']*Infection-Sigma()*Midpar1.loc[0,c+h+w+'_E1']-imigration_rate*Midpar1.loc[0,c+h+w+'_E1']+qprob*emigration_num_E)\n",
    "                Midpar2.loc[0,c+h+w+'_I2']=Midpar1.loc[0,c+h+w+'_I1']+dt*(Sigma()*Midpar1.loc[0,c+h+w+'_E1']-(Recover+Dead)*Midpar1.loc[0,c+h+w+'_I1']-imigration_rate*Midpar1.loc[0,c+h+w+'_I1']+qprob*emigration_num_I)\n",
    "                Midpar2.loc[0,c+h+w+'_R2']=Midpar1.loc[0,c+h+w+'_R1']+dt*(Recover*Midpar1.loc[0,c+h+w+'_I1']-imigration_rate*Midpar1.loc[0,c+h+w+'_R1']+qprob*emigration_num_R)\n",
    "                Midpar2.loc[0,c+h+w+'_D2']=Midpar1.loc[0,c+h+w+'_D1']+dt*Dead*Midpar1.loc[0,c+h+w+'_I1']\n",
    "                Midpar2.loc[0,c+h+w+'_H2']=Midpar1.loc[0,c+h+w+'_H1']+dt*Sigma()*Midpar1.loc[0,c+h+w+'_E1']\n",
    "                Midpar1.loc[0,c+h+w+'_S1']=Midpar2.loc[0,c+h+w+'_S2']\n",
    "                Midpar1.loc[0,c+h+w+'_E1']=Midpar2.loc[0,c+h+w+'_E2']\n",
    "                Midpar1.loc[0,c+h+w+'_I1']=Midpar2.loc[0,c+h+w+'_I2']\n",
    "                Midpar1.loc[0,c+h+w+'_R1']=Midpar2.loc[0,c+h+w+'_R2']\n",
    "                Midpar1.loc[0,c+h+w+'_D1']=Midpar2.loc[0,c+h+w+'_D2']\n",
    "                Midpar1.loc[0,c+h+w+'_H1']=Midpar2.loc[0,c+h+w+'_H2']\n",
    "                Solution.loc[i,c+h+w+'_S']=Midpar2.loc[0,c+h+w+'_S2']\n",
    "                Solution.loc[i,c+h+w+'_E']=Midpar2.loc[0,c+h+w+'_E2']\n",
    "                Solution.loc[i,c+h+w+'_I']=Midpar2.loc[0,c+h+w+'_I2']\n",
    "                Solution.loc[i,c+h+w+'_R']=Midpar2.loc[0,c+h+w+'_R2']\n",
    "                Solution.loc[i,c+h+w+'_D']=Midpar2.loc[0,c+h+w+'_D2']\n",
    "                Solution.loc[i,c+h+w+'_H']=Midpar2.loc[0,c+h+w+'_H2']\n",
    "#输出数据\n",
    "#S^{c,h}\n",
    "#每个城市各区感染人口\n",
    "INDEX_Distinct=[]\n",
    "for  c in Citys:\n",
    "    for h in Patches[c]:\n",
    "        INDEX_Distinct.append(c+h)\n",
    "COLUMNS_Distinct=[] #每列\n",
    "for i in INDEX_Distinct:\n",
    "    A=[i+'_S',i+'_E',i+'_I',i+'_R',i+'_D',i+'_H']\n",
    "    COLUMNS_Distinct.extend(A)\n",
    "Solution_Distinct=pd.DataFrame(columns=COLUMNS_Distinct)\n",
    "for i in range(int(Tspan/dt)):\n",
    "    for c in Citys:\n",
    "        for h in Patches[c]:\n",
    "            S_ch=0\n",
    "            E_ch=0\n",
    "            I_ch=0\n",
    "            R_ch=0\n",
    "            D_ch=0\n",
    "            H_ch=0\n",
    "            for w in Patches[c]:\n",
    "                S_ch=S_ch+Solution.loc[i,c+h+w+'_S']\n",
    "                E_ch=E_ch+Solution.loc[i,c+h+w+'_E']\n",
    "                I_ch=I_ch+Solution.loc[i,c+h+w+'_I']\n",
    "                R_ch=R_ch+Solution.loc[i,c+h+w+'_R']\n",
    "                D_ch=D_ch+Solution.loc[i,c+h+w+'_D']\n",
    "                H_ch=H_ch+Solution.loc[i,c+h+w+'_H']\n",
    "            Solution_Distinct.loc[i,c+h+'_S']=S_ch\n",
    "            Solution_Distinct.loc[i,c+h+'_E']=E_ch\n",
    "            Solution_Distinct.loc[i,c+h+'_I']=I_ch\n",
    "            Solution_Distinct.loc[i,c+h+'_R']=R_ch\n",
    "            Solution_Distinct.loc[i,c+h+'_D']=D_ch\n",
    "            Solution_Distinct.loc[i,c+h+'_H']=H_ch\n",
    "#各市感染人口\n",
    "\n",
    "COLUMNS_City=[] #每列\n",
    "for i in Citys:\n",
    "    A=[i+'_S',i+'_E',i+'_I',i+'_R',i+'_D',i+'_H']\n",
    "    COLUMNS_City.extend(A)\n",
    "Solution_City=pd.DataFrame(columns=COLUMNS_City)           \n",
    "for i in range(int(Tspan/dt)):\n",
    "    for c in Citys:\n",
    "        S_c=0\n",
    "        E_c=0\n",
    "        I_c=0\n",
    "        R_c=0\n",
    "        D_c=0\n",
    "        H_c=0        \n",
    "        for h in Patches[c]:\n",
    "            S_c=S_c+Solution_Distinct.loc[i,c+h+'_S']\n",
    "            E_c=E_c+Solution_Distinct.loc[i,c+h+'_E']\n",
    "            I_c=I_c+Solution_Distinct.loc[i,c+h+'_I']\n",
    "            R_c=R_c+Solution_Distinct.loc[i,c+h+'_R']\n",
    "            D_c=D_c+Solution_Distinct.loc[i,c+h+'_D']\n",
    "            H_c=H_c+Solution_Distinct.loc[i,c+h+'_H']\n",
    "        Solution_City.loc[i,c+'_S']=S_c\n",
    "        Solution_City.loc[i,c+'_E']=E_c\n",
    "        Solution_City.loc[i,c+'_I']=I_c\n",
    "        Solution_City.loc[i,c+'_R']=R_c\n",
    "        Solution_City.loc[i,c+'_D']=D_c\n",
    "        Solution_City.loc[i,c+'_H']=H_c\n",
    "            \n",
    "print(Solution_City)\n",
    "end=time.time()\n",
    "print('Running time: %s Seconds'%(end-start))\n",
    "\n",
    "#保存数据\n",
    "Solution_City.to_csv('Down95City202003302.csv',index=0) #市\n",
    "Solution_Distinct.to_csv('Down95Distinct202003302.csv',index=0)#区\n",
    "Solution.to_csv('Down95Detail202003302.csv',index=0)#人群子类\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13\n",
      "<class 'list'>\n",
      "Help on class Map in module pyecharts.charts.basic_charts.map:\n",
      "\n",
      "class Map(pyecharts.charts.chart.Chart, MapMixin)\n",
      " |  Map(init_opts: Union[pyecharts.options.global_options.InitOpts, dict] = <pyecharts.options.global_options.InitOpts object at 0x00000219185FE348>)\n",
      " |  \n",
      " |  `Base` is the root class for all graphical class, it provides\n",
      " |  part of the initialization parameters and common methods\n",
      " |  \n",
      " |  Method resolution order:\n",
      " |      Map\n",
      " |      pyecharts.charts.chart.Chart\n",
      " |      pyecharts.charts.base.Base\n",
      " |      pyecharts.charts.mixins.ChartMixin\n",
      " |      MapMixin\n",
      " |      builtins.object\n",
      " |  \n",
      " |  Methods inherited from pyecharts.charts.chart.Chart:\n",
      " |  \n",
      " |  __init__(self, init_opts: Union[pyecharts.options.global_options.InitOpts, dict] = <pyecharts.options.global_options.InitOpts object at 0x00000219185FE348>)\n",
      " |      Initialize self.  See help(type(self)) for accurate signature.\n",
      " |  \n",
      " |  set_colors(self, colors: Sequence[str])\n",
      " |  \n",
      " |  set_global_opts(self, title_opts: Union[pyecharts.options.global_options.TitleOpts, dict] = <pyecharts.options.global_options.TitleOpts object at 0x0000021918608088>, legend_opts: Union[pyecharts.options.global_options.LegendOpts, dict] = <pyecharts.options.global_options.LegendOpts object at 0x0000021918608148>, tooltip_opts: Union[pyecharts.options.global_options.TooltipOpts, dict, NoneType] = None, toolbox_opts: Union[pyecharts.options.global_options.ToolboxOpts, dict] = None, brush_opts: Union[pyecharts.options.global_options.BrushOpts, dict, NoneType] = None, xaxis_opts: Union[pyecharts.options.global_options.AxisOpts, dict, NoneType] = None, yaxis_opts: Union[pyecharts.options.global_options.AxisOpts, dict, NoneType] = None, visualmap_opts: Union[pyecharts.options.global_options.VisualMapOpts, dict, Sequence[Union[pyecharts.options.global_options.VisualMapOpts, dict]], NoneType] = None, datazoom_opts: Union[pyecharts.options.global_options.DataZoomOpts, dict, Sequence[Union[pyecharts.options.global_options.DataZoomOpts, dict]], NoneType] = None, graphic_opts: Union[pyecharts.options.charts_options.BaseGraphic, dict, Sequence[Union[pyecharts.options.charts_options.BaseGraphic, dict]], NoneType] = None, axispointer_opts: Union[pyecharts.options.global_options.AxisPointerOpts, dict, NoneType] = None)\n",
      " |  \n",
      " |  set_series_opts(self, label_opts: Union[pyecharts.options.series_options.LabelOpts, dict, NoneType] = None, linestyle_opts: Union[pyecharts.options.series_options.LineStyleOpts, dict, NoneType] = None, splitline_opts: Union[pyecharts.options.series_options.SplitLineOpts, dict, NoneType] = None, areastyle_opts: Union[pyecharts.options.series_options.AreaStyleOpts, dict, NoneType] = None, axisline_opts: Union[pyecharts.options.global_options.AxisLineOpts, dict, NoneType] = None, markpoint_opts: Union[pyecharts.options.series_options.MarkPointOpts, dict, NoneType] = None, markline_opts: Union[pyecharts.options.series_options.MarkLineOpts, dict, NoneType] = None, markarea_opts: Union[pyecharts.options.series_options.MarkAreaOpts, dict, NoneType] = None, effect_opts: Union[pyecharts.options.series_options.EffectOpts, dict, NoneType] = <pyecharts.options.series_options.EffectOpts object at 0x00000219185FE3C8>, tooltip_opts: Union[pyecharts.options.global_options.TooltipOpts, dict, NoneType] = None, itemstyle_opts: Union[pyecharts.options.series_options.ItemStyleOpts, dict, NoneType] = None, **kwargs)\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Methods inherited from pyecharts.charts.base.Base:\n",
      " |  \n",
      " |  dump_options(self) -> str\n",
      " |  \n",
      " |  dump_options_with_quotes(self) -> str\n",
      " |  \n",
      " |  get_options(self) -> dict\n",
      " |  \n",
      " |  render(self, path: str = 'render.html', template_name: str = 'simple_chart.html', env: Union[jinja2.environment.Environment, NoneType] = None, **kwargs) -> str\n",
      " |  \n",
      " |  render_embed(self, template_name: str = 'simple_chart.html', env: Union[jinja2.environment.Environment, NoneType] = None, **kwargs) -> str\n",
      " |  \n",
      " |  render_notebook(self)\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Methods inherited from pyecharts.charts.mixins.ChartMixin:\n",
      " |  \n",
      " |  add_js_funcs(self, *fns)\n",
      " |  \n",
      " |  load_javascript(self)\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Data descriptors inherited from pyecharts.charts.mixins.ChartMixin:\n",
      " |  \n",
      " |  __dict__\n",
      " |      dictionary for instance variables (if defined)\n",
      " |  \n",
      " |  __weakref__\n",
      " |      list of weak references to the object (if defined)\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Methods inherited from MapMixin:\n",
      " |  \n",
      " |  add(self, series_name: str, data_pair: Sequence, maptype: str = 'china', *, is_selected: bool = True, is_roam: bool = True, center: Union[Sequence, NoneType] = None, zoom: Union[int, float, NoneType] = 1, name_map: Union[dict, NoneType] = None, symbol: Union[str, NoneType] = None, is_map_symbol_show: bool = True, label_opts: Union[pyecharts.options.series_options.LabelOpts, dict, NoneType] = <pyecharts.options.series_options.LabelOpts object at 0x000002191861D408>, tooltip_opts: Union[pyecharts.options.global_options.TooltipOpts, dict, NoneType] = None, itemstyle_opts: Union[pyecharts.options.series_options.ItemStyleOpts, dict, NoneType] = None, emphasis_label_opts: Union[pyecharts.options.series_options.LabelOpts, dict, NoneType] = None, emphasis_itemstyle_opts: Union[pyecharts.options.series_options.ItemStyleOpts, dict, NoneType] = None)\n",
      "\n",
      "None\n"
     ]
    }
   ],
   "source": [
    "from pyecharts.charts import Map\n",
    "from pyecharts import options as opts\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "def map_wuhan(data,Cmax):\n",
    "    c = (\n",
    "        Map()\n",
    "        .add(\"武汉-Map\", data, \"武汉\")\n",
    "        .set_global_opts(\n",
    "            title_opts=opts.TitleOpts(title=\"Map-武汉地图\"),\n",
    "            visualmap_opts=opts.VisualMapOpts(max_=Cmax),\n",
    "        )\n",
    "    )\n",
    "    return c\n",
    "def map_beijing(data,Cmax):\n",
    "    c = (\n",
    "        Map()\n",
    "        .add(\"北京-Map\", data, \"北京\")\n",
    "        .set_global_opts(\n",
    "            title_opts=opts.TitleOpts(title=\"Map-北京地图\"),\n",
    "            visualmap_opts=opts.VisualMapOpts(max_=Cmax),\n",
    "        )\n",
    "    )\n",
    "    return c \n",
    "\n",
    "Wuhanqu=['江岸区','江汉区','硚口区','汉阳区','武昌区','青山区',\n",
    "         '洪山区','蔡甸区','江夏区','黄陂区','新洲区','东西湖区','汉南区']\n",
    "Beijingqu=['密云区','延庆区','朝阳区','丰台区','石景山区','海淀区','门头沟区', \n",
    "           '房山区','通州区','顺义区','昌平区','大兴区','怀柔区','平谷区','东城区','西城区']\n",
    "values3= np.array([4300, 7290, 7239, 3579, 8224, 2958,4999,2588,1974,1661,1775,1018,1517])\n",
    "\n",
    "values_renkou=np.array([96.24,72.96,86.85,65.27,127.63,52.88,163.75,73.50,91.37,98.83,90.21,56.25,13.55])\n",
    "print(len(values_renkou))\n",
    "values_severity=values3/values_renkou\n",
    "values2=[8,9,20,3, 5, 7, 8, 2, 4,10,11,25,26,40,1,5]\n",
    "\n",
    "data_H = []\n",
    "for index in range(len(Wuhanqu)):\n",
    "    city_info = [Wuhanqu[index], values3[index]]\n",
    "    data_H.append(city_info)\n",
    "\n",
    "data_wuhanrenkou = []\n",
    "for index in range(len(Wuhanqu)):\n",
    "    city_info = [Wuhanqu[index], values_renkou[index]]\n",
    "    data_wuhanrenkou.append(city_info)\n",
    "data_sevirity = []\n",
    "for index in range(len(Wuhanqu)):\n",
    "    city_info = [Wuhanqu[index], values_severity[index]]\n",
    "    data_sevirity.append(city_info)\n",
    "    \n",
    "data_B = []\n",
    "for index in range(len(Beijingqu)):\n",
    "    city_info = [Beijingqu[index], values2[index]]\n",
    "    data_B.append(city_info)\n",
    "\n",
    "print(type(data_B))\n",
    "\n",
    "\n",
    "a = map_wuhan(data_wuhanrenkou,max(values_renkou))\n",
    "b = map_beijing(data_B,50)\n",
    "c =  map_wuhan(data_sevirity,max(values_severity))\n",
    "a.render(path=\"./mapresults/武汉人口地图.html\")\n",
    "b.render(path=\"./mapresults/beijing地图.html\")\n",
    "c.render(path=\"./mapresults/武汉疫情严重程度地图.html\")\n",
    "print(help(Map))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "16\n"
     ]
    }
   ],
   "source": [
    "Citys=['Beijing','Wuhan']\n",
    "Wuhanqu=['江岸区','江汉区','硚口区','汉阳区','武昌区','青山区',\n",
    "         '洪山区','蔡甸区','江夏区','黄陂区','新洲区','东西湖区','汉南区']\n",
    "Beijingqu=['密云区','延庆区','朝阳区','丰台区','石景山区','海淀区','门头沟区', \n",
    "           '房山区','通州区','顺义区','昌平区','大兴区','怀柔区','平谷区','东城区','西城区']\n",
    "A={'Beijing':Beijingqu,'Wuhan':Wuhanqu}\n",
    "print(len(A['Beijing']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "print('Beijing' in 'Beijings')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<generator object <genexpr> at 0x0000026C619DAC78>\n"
     ]
    }
   ],
   "source": [
    "Patches={'Beijing':Beijingqu,'Wuhan':Wuhanqu}\n",
    "Wuhanqu=['江岸区','江汉区','硚口区','汉阳区','武昌区','青山区',\n",
    "         '洪山区','蔡甸区','江夏区','黄陂区','新洲区','东西湖区','汉南区']\n",
    "Beijingqu=['密云区','延庆区','朝阳区','丰台区','石景山区','海淀区','门头沟区', \n",
    "           '房山区','通州区','顺义区','昌平区','大兴区','怀柔区','平谷区','东城区','西城区']\n",
    "print(len(Patches[i]) for i in ['Beijing','Wuhan'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n"
     ]
    }
   ],
   "source": [
    "a=1+-1\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11.756301359092745\n",
      "16678.076923076926\n"
     ]
    }
   ],
   "source": [
    "print(216815/(13*(908.35+510.30)))\n",
    "print(11.756301359092745*(908.35+510.30))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "print(1>2 or 2>1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
